#####################################################################
# Table 6
#####################################################################

#####################################################################
# Preliminaries
#####################################################################

rm(list = ls())

library(xtable)
library(scinference)
library(limSolve)
library(writexl)

#####################################################################
# Functions
#####################################################################

source("code/auxiliary scripts/common_functions.R")

#####################################################################
# Calibration DGPs
#####################################################################

source("code/auxiliary scripts/calibration_dgps.R")

#####################################################################
# Simulations for large T0
#####################################################################

# Setup

set.seed(12345)

nreps_ttest <- 10000

# Simulations 

print("LaTeX output corresponding to Table 6, printed separately for each DGP")

results_overall <- matrix(NA,9,9)

for (DGP in 1:9){
  
  begin <- Sys.time()
  
  results_K4_cov <- results_K4_leng <- results_K4_bias <- matrix(NA,nreps_ttest,3)
  results_K6_cov <- results_K6_leng <- results_K6_bias <- matrix(NA,nreps_ttest,3)
  results_K8_cov <- results_K8_leng <- results_K8_bias <- matrix(NA,nreps_ttest,3)

  for (r in 1:nreps_ttest){
    res_K4_temp    <- sim_one_sample(DGP,T0=150,T1=T1_app,J=J_app,K=4,Lambda,rho_u,var_u,var_factors,rho_vec,var_epsl_vec,w0_sc,tau_alternative=0,li=FALSE,sdid=FALSE,permtest = FALSE)
    res_K6_temp    <- sim_one_sample(DGP,T0=150,T1=T1_app,J=J_app,K=6,Lambda,rho_u,var_u,var_factors,rho_vec,var_epsl_vec,w0_sc,tau_alternative=0,li=FALSE,sdid=FALSE,permtest = FALSE)
    res_K8_temp    <- sim_one_sample(DGP,T0=150,T1=T1_app,J=J_app,K=8,Lambda,rho_u,var_u,var_factors,rho_vec,var_epsl_vec,w0_sc,tau_alternative=0,li=FALSE,sdid=FALSE,permtest = FALSE)

    results_K4_bias[r,]  <- res_K4_temp$bias_all
    results_K4_leng[r,]  <- res_K4_temp$leng_all
    results_K4_cov[r,]   <- res_K4_temp$cov_all
    
    results_K6_bias[r,]  <- res_K6_temp$bias_all
    results_K6_leng[r,]  <- res_K6_temp$leng_all
    results_K6_cov[r,]   <- res_K6_temp$cov_all
    
    results_K8_bias[r,]  <- res_K8_temp$bias_all
    results_K8_leng[r,]  <- res_K8_temp$leng_all
    results_K8_cov[r,]   <- res_K8_temp$cov_all
    
  }
  
  cf_K4_bias <- colMeans(results_K4_bias)
  cf_K6_bias <- colMeans(results_K6_bias)
  cf_K8_bias <- colMeans(results_K8_bias)

  cf_K4_cov <- colMeans(results_K4_cov)
  cf_K6_cov <- colMeans(results_K6_cov)
  cf_K8_cov <- colMeans(results_K8_cov)  

  cf_K4_leng <- colMeans(results_K4_leng)
  cf_K6_leng <- colMeans(results_K6_leng)
  cf_K8_leng <- colMeans(results_K8_leng)  

  print(DGP)
  print(Sys.time()-begin)
  results_table <- c(10*cf_K4_bias[1],10*cf_K6_bias[1],10*cf_K8_bias[1],cf_K4_cov[1],cf_K6_cov[1],cf_K8_cov[1],cf_K4_leng[1],cf_K6_leng[1],cf_K8_leng[1])
  print(xtable(rbind(results_table),digits=c(rep(2,10))),include.rownames = F)
  
  results_overall[DGP,] <- sprintf("%.2f",results_table)
  
}

# Save table in .xlsx format

data_results_table6 <- as.data.frame(results_overall)
colnames(data_results_table6) <- NULL
write_xlsx(data_results_table6,"tables/table6.xlsx")

